!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!> \file
!! Main program for the CG method, Stokes problem.
!!
!! \n
!!
!! The input file is composed by one <em>namelist</em> named \c input
!! as described in \c mod_cg_ns_setup. The name of the input file can
!! be specified as command argument, otherwise the default \c
!! input_file_name_def is used.
!!
!! Other relevant modules are
!! <ul>
!!  <li> \c mod_cg_ns_linsystem: representation of the linear system
!! </ul>
!<
program cg_stokes_main

 use mod_utils, only: &
   t_realtime, my_second

 use mod_messages, only: &
   mod_messages_constructor, &
   mod_messages_destructor,  &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   wp

 use mod_fu_manager, only: &
   new_file_unit

 use mod_octave_io, only: &
   write_octave

 use mod_sparse, only: &
   ! sparse types
   t_col,       &
   t_tri,       &
   ! overloaded operators
   matmul,      &
   clear

 use mod_octave_io_sparse, only: &
   write_octave

 use mod_linsolver, only: &
   factor, solve, clean

 use mod_output_control, only: &
   elapsed_format, &
   base_name

 use mod_grid, only: &
   t_grid,  &
   grid

 use mod_cgdofs, only: &
   t_cgdofs, new_cgdofs, &
   write_octave

 use mod_base, only: &
   t_base,       &
   write_octave

 use mod_testcases, only: &
   test_name,   &
   test_description,&
   coeff_dir,   &
   coeff_u,     &
   coeff_gradu, &
   coeff_p,     &
   coeff_gradp

 use mod_cg_ns_setup, only: &
   mod_cg_ns_setup_constructor, &
   mod_cg_ns_setup_destructor,  &
   ! miscellaneous
   max_char_len,   &
   ! FE space
   pressure_stabilization, add_div_lagmultiplier, &
   ku, kp,         &
   p_stab_coeff,   &
   uref_nodes, pref_nodes, &
   ubase, pbase,   &
   udofs, pdofs,   &
   ! linear solver
   linsolver,      &
   umfpack_print_level, &
   gmres_abstol,   &
   precon_ndiags,  &
   gmres_toll,     &
   preconditioner, &
   ! error computation
   compute_error_norms, &
   err_extra_deg,  &
   ! IO
   write_sys

 use mod_cg_ns_locmat, only: &
   t_bcs_error

 use mod_cg_ns_linsystem, only: &
   mod_cg_ns_linsystem_constructor, &
   mod_cg_ns_linsystem_destructor,  &
   ns_mat, ns_partmat, ns_packsol

 use mod_error_norms, only: &
   mod_error_norms_constructor, &
   mod_error_norms_destructor,  &
   velocity_err, pressure_err

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------
 
 ! Parameters
 character(len=*), parameter :: input_file_name_def = 'cg-stokes.in'
 character(len=*), parameter :: this_prog_name = 'cg-stokes'

 ! Local matrixes
 logical :: static_condensation
 real(wp), allocatable :: a22ia21(:,:,:), a22ib2t(:,:,:), a22if2(:,:)
 type(t_bcs_error) :: b_err

 ! Global matrix
 real(wp), allocatable :: fff(:), fff1(:), fff2(:), rhs(:), &
  uuu(:), uuu1(:), uuu2(:)
 type(t_col) :: mmm, mmm11, mmm12, mmm21, mmm22

 ! indexes
 integer :: ie, i

 ! bcs variables
 integer, allocatable :: dofs_nat(:), dofs_dir(:)

 ! error computations
 real(wp), allocatable :: uuu_dg(:,:), ppp_dg(:)
 real(wp) :: e_u_l2, e_u_h1, e_u_div, e_p_l2, e_p_h1

 ! IO variables
 character(len=*), parameter :: out_file_results_suff = '-results.octxt'
 character(len=10000+len(out_file_results_suff)) :: out_file_results_name

 ! Auxiliary variables
 integer :: fu, ierr
 real(t_realtime) :: t00, t0, t1
 character(len=max_char_len) :: message

 !-----------------------------------------------------------------------
 ! Initializations and stratup
 t00 = my_second()

 call mod_messages_constructor()

 call mod_cg_ns_setup_constructor(input_file_name_def)
 !-----------------------------------------------------------------------

 !-----------------------------------------------------------------------
 ! Assemble matrix and right hand side
 call mod_cg_ns_linsystem_constructor(ubase,pbase,p_stab_coeff)

 t0 = my_second()
 call ns_mat(mmm,fff,static_condensation,a22ia21,a22ib2t,a22if2, &
             uref_nodes,ubase,pbase,udofs,pdofs,                 &
             add_div_lagmultiplier,pressure_stabilization,b_err)

 if(b_err%lerr) call error(this_prog_name,'',b_err%message)
 t1 = my_second()
 write(message,elapsed_format) &
   'Completed assembling of the linear system: elapsed time ',t1-t0,'s.'
 call info(this_prog_name,'',message)
 !-----------------------------------------------------------------------

 !-----------------------------------------------------------------------
 ! Partition mmm and the rhs according to the Dirichlet bcs.
 t0 = my_second()

 call ns_partmat(mmm11,mmm12,mmm21,mmm22,fff1,uuu2,mmm,fff,     &
   udofs,pdofs,add_div_lagmultiplier,dofs_nat,dofs_dir,coeff_dir)

 t1 = my_second()
 write(message,elapsed_format) &
   'Completed partitioning of the linear system: elapsed time ',t1-t0,'s.'
 call info(this_prog_name,'',message)
 !-----------------------------------------------------------------------

 !-----------------------------------------------------------------------
 ! Solution of the linear system
 t0 = my_second()

 allocate( uuu1(mmm11%n) , rhs(mmm11%n) )

 rhs = fff1 - matmul(mmm12,uuu2)

 call factor(mmm11)
 call solve(mmm11,uuu1,rhs)
 call clean()

 ! boundary terms
 allocate( fff2(mmm22%n) )
 fff2 = matmul(mmm21,uuu1) + matmul(mmm22,uuu2)
 fff(dofs_dir+1) = fff2

 deallocate( rhs , fff1 , fff2)
 call clear( mmm11 )
 call clear( mmm12 )
 call clear( mmm21 )
 call clear( mmm22 )

 t1 = my_second()
 write(message,elapsed_format) &
   'Completed solution of the linear system: elapsed time ',t1-t0,'s.'
 call info(this_prog_name,'',message)
 !-----------------------------------------------------------------------

 !-----------------------------------------------------------------------
 ! Postprocessings
 t1 = my_second()

 call ns_packsol(uuu,uref_nodes,ubase,pbase,udofs,pdofs, &
                 dofs_nat,dofs_dir,add_div_lagmultiplier,&
                 uuu1,uuu2,a22if2,a22ia21,a22ib2t)

 deallocate( uuu1 , uuu2 , a22if2 , a22ia21 , a22ib2t )

 t1 = my_second()
 write(message,elapsed_format) &
   'Completed post-processings: elapsed time ',t1-t0,'s.'
 call info(this_prog_name,'',message)
 !-----------------------------------------------------------------------

 !-----------------------------------------------------------------------
 ! Error computation
 t0 = my_second()

 ! In order to re-use as much as possible the error computation
 ! subroutines written for the DG case, we reshape the numerical
 ! solution in terms of "discontinuous" functions, by formally
 ! assigning different degrees of freedom to each elements. Of course,
 ! these additional degrees of freedom are indeed copies of the CG
 ! dofs.
 allocate(uuu_dg(grid%ne*ubase%pk,grid%d),ppp_dg(grid%ne*pbase%pk))

 do ie=1,grid%ne
   do i=1,grid%d
     uuu_dg( (ie-1)*ubase%pk+1 : ie*ubase%pk , i ) = &
       uuu( (udofs%dofs(:,ie)-1)*grid%d + i )
   enddo
   ppp_dg( (ie-1)*pbase%pk+1 : ie*pbase%pk ) = &
     uuu( grid%d*udofs%ndofs + pdofs%dofs(:,ie) )
 enddo
 errnorms: if(compute_error_norms) then
   call mod_error_norms_constructor()

   call velocity_err(ubase,uuu_dg,coeff_u,coeff_gradu,err_extra_deg, &
                     e_u_l2, e_u_h1, e_u_div)
   call pressure_err(pbase,ppp_dg,coeff_p,coeff_gradp,err_extra_deg, &
                     e_p_l2, e_p_h1)

   call mod_error_norms_destructor()
 endif errnorms

 t1 = my_second()
 write(message,elapsed_format) &
   'Completed error computations: elapsed time ',t1-t0,'s.'
 call info(this_prog_name,'',message)
 !-----------------------------------------------------------------------

 !-----------------------------------------------------------------------
 ! Write the results
 call new_file_unit(fu,ierr)
 out_file_results_name = trim(base_name)//out_file_results_suff
 open(fu,file=trim(out_file_results_name), &
      status='replace',action='write',form='formatted',iostat=ierr)
   ! preamble
   call write_octave(test_name       ,'test_name'       ,fu)
   call write_octave(test_description,'test_description',fu)
   ! data
   if(write_sys) then
     call write_octave(dofs_nat,'c','dofs_nat',fu)
     call write_octave(dofs_dir,'c','dofs_dir',fu)
     call write_octave(mmm,'mmm',fu)
     call write_octave(fff,'c','fff',fu)
   endif
   call write_octave(uuu,'c','uuu',fu)
   call write_octave(uuu_dg,'uuu_dg',fu)
   call write_octave(ppp_dg,'c','ppp_dg',fu)
   if(compute_error_norms) then
     call write_octave(e_u_l2 ,'e_u_l2' ,fu)
     call write_octave(e_u_h1 ,'e_u_h1' ,fu)
     call write_octave(e_u_div,'e_u_div',fu)
     call write_octave(e_p_l2 ,'e_p_l2' ,fu)
     call write_octave(e_p_h1 ,'e_p_h1' ,fu)
   endif
 close(fu,iostat=ierr)
 !-----------------------------------------------------------------------

 !-----------------------------------------------------------------------
 ! Cleanup
 deallocate( dofs_nat , dofs_dir , fff , uuu )
 deallocate( uuu_dg , ppp_dg )
 call clear( mmm )

 call mod_cg_ns_linsystem_destructor()

 call mod_cg_ns_setup_destructor()

 t1 = my_second()
 write(message,elapsed_format) &
   'Done: elapsed time ',t1-t00,'s.'
 call info(this_prog_name,'',message)

 call mod_messages_destructor()
 !-----------------------------------------------------------------------

end program cg_stokes_main

